{
 "metadata": {
  "name": "",
  "signature": "sha256:5f4abf83f0e42e9967bfabd0e023e8f806282c77b69a3bf7926de11defc6314e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Adding a new optimization problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this Tutorial we will learn how to code simple optimization problems (continuous, single objective, unconstrained), so that PyGMO can then apply all of its algorithmic power to solve it. In a nutshell, we will write a class deriving from **PyGMO.problem.base** and reimplement some of its \u2018virtual\u2019 methods."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Single objective problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us start with defining one of the classic textbook examples of an optimization problem."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from PyGMO.problem import base\n",
      "\n",
      "\n",
      "class my_problem(base):\n",
      "    \"\"\"\n",
      "    De Jong (sphere) function implemented purely in Python.\n",
      "\n",
      "    USAGE: my_problem(dim=10)\n",
      "\n",
      "    * dim problem dimension\n",
      "    \"\"\"\n",
      "\n",
      "    def __init__(self, dim=10):\n",
      "        # First we call the constructor of the base class telling PyGMO\n",
      "        # what kind of problem to expect ('dim' dimensions, 1 objective, 0 contraints etc.)\n",
      "        super(my_problem, self).__init__(dim)\n",
      "\n",
      "        # We set the problem bounds (in this case equal for all components)\n",
      "        self.set_bounds(-5.12, 5.12)\n",
      "\n",
      "    # Reimplement the virtual method that defines the objective function.\n",
      "    def _objfun_impl(self, x):\n",
      "\n",
      "        # Compute the sphere function\n",
      "        f = sum([x[i] ** 2 for i in range(self.dimension)])\n",
      "\n",
      "        # Note that we return a tuple with one element only. In PyGMO the objective functions\n",
      "        # return tuples so that multi-objective optimization is also possible.\n",
      "        return (f, )\n",
      "\n",
      "    # Finally we also reimplement a virtual method that adds some output to the __repr__ method\n",
      "    def human_readable_extra(self):\n",
      "        return \"\\n\\t Problem dimension: \" + str(self.__dim)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that by default PyGMO will assume one wants to minimize the objective function. In the second part of this tutorial we will also see how it is possible to change this default behaviour.\n",
      "\n",
      "To solve our problem we will use Artificial Bee Colony algorithm with 20 individuals."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from PyGMO import algorithm, island\n",
      "\n",
      "prob = my_problem(dim=10)  # Create a 10-dimensional problem\n",
      "algo = algorithm.bee_colony(gen=500)  # 500 generations of bee_colony algorithm\n",
      "isl = island(algo, prob, 20)  # Instantiate population with 20 individuals\n",
      "isl.evolve(1)  # Evolve the island once\n",
      "isl.join()\n",
      "print(isl.population.champion.f)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And we are done! Objective value in the order of $10^{-25}$, no big deal for a sphere problem.\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Single objective maximization problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let\u2019s consider now a maximization problem. To solve such a problem, two possibilities are available to the PaGMO/PyGMO user. The first one is to code the original problem as a minimization problem by premultiplying the objective function by $-1$ (a technique wich is often used and requires no particular effort). If such a method is used, the final fitness value obtained with PyGMO has to be multiplied by $-1$ to get back to the correct value.\n",
      "\n",
      "A second method, more elegant and most of all serving the purpose to show the use of another virtual method which can be reimplemented in python objects deriving from base, is to override the function that compares two fitness vectors. This function is used by all pagmo algorithms to compare performances of individuals. By default, this function compares the fitness $f_1$ to a fitness $f_2$ and returns true if $f_1$ dominates $f_2$ (which is single objective optimization correspond to minimization). Let us see how..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class my_problem_max(base):\n",
      "    \"\"\"\n",
      "    Analytical function to maximize.\n",
      "\n",
      "    USAGE: my_problem_max()\n",
      "    \"\"\"\n",
      "    \n",
      "    def __init__(self):\n",
      "        super(my_problem_max, self).__init__(2)\n",
      "        self.set_bounds(-10, 10)\n",
      "        \n",
      "        # We provide a list of the best known solutions to the problem\n",
      "        self.best_x = [[1.0, -1.0], ]\n",
      "\n",
      "    # Reimplement the virtual method that defines the objective function\n",
      "    def _objfun_impl(self, x):\n",
      "        f = -(1.0 - x[0]) ** 2 - 100 * (-x[0] ** 2 - x[1]) ** 2 - 1.0\n",
      "        return (f, )\n",
      "\n",
      "    # Reimplement the virtual method that compares fitnesses\n",
      "    def _compare_fitness_impl(self, f1, f2):\n",
      "        return f1[0] > f2[0]\n",
      "\n",
      "    # Add some output to __repr__\n",
      "    def human_readable_extra(self):\n",
      "        return \"\\n\\tMaximization problem\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Additionally in the constructor we provide a list of all known global minima (we will use those later for testing).\n",
      "The list of corresponding objective function values will be then computed and accessible through **best_f** of the problem's instance.\n",
      "\n",
      "As before, we use our favorite optimization algorithm:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from math import sqrt\n",
      "\n",
      "prob = my_problem_max()\n",
      "algo = algorithm.de(gen=20)\n",
      "isl = island(algo, prob, 20)\n",
      "isl.evolve(10)\n",
      "isl.join()\n",
      "\n",
      "print(\"Best individual:\")\n",
      "print(isl.population.champion)\n",
      "\n",
      "print(\"Comparison of the best found fitness with the best known fitness:\")\n",
      "for best_fitness in prob.best_f:\n",
      "    print(best_fitness[0] - isl.population.champion.f[0])\n",
      "\n",
      "print(\"L2 distance to the best decision vector:\")\n",
      "for best_decision in prob.best_x:\n",
      "    l2_norm = 0\n",
      "    for n in range(0, len(best_decision)):\n",
      "        l2_norm +=  (best_decision[n] - isl.population.champion.x[n]) ** 2\n",
      "    l2_norm = sqrt(l2_norm)\n",
      "    print(l2_norm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note here that we used the **best_f** and **best_x** methods which return the best known fitness and decision vectors. The **best_f** vector is automatically available as we defined **best_x** in the problem. With these vectors, we can have an idea of the optimizer performances. The result of this optimization should be in order of $10^{-11}$ for the comparison with the best fitness and $10^{-6}$ for the distance to the best decision vector."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Multi-objective problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As hinted before, users can also define their own multi-objective problem.\n",
      "In that case we need to overload the the base constructor with third argument stating the desired objective function dimension and return a tuple or a list with more than one element in the objective function implementation (both dimensions must agree)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class my_mo_problem(base):\n",
      "    \"\"\"\n",
      "    A multi-objective problem.\n",
      "    (This is actually a Python implementation of 2-dimensional ZDT-1 problem)\n",
      "\n",
      "    USAGE: my_mo_problem()\n",
      "    \"\"\"\n",
      "    \n",
      "    def __init__(self, dim=2):\n",
      "        # We call the base constructor as 'dim' dimensional problem, with 0 integer parts and 2 objectives.\n",
      "        super(my_mo_problem,self).__init__(dim, 0, 2)\n",
      "        self.set_bounds(0.0, 1.0)\n",
      "\n",
      "    # Reimplement the virtual method that defines the objective function\n",
      "    def _objfun_impl(self, x):\n",
      "        f0 = x[0]\n",
      "        g = 1.0 + 4.5 * x[1]\n",
      "        f1 = g * (1.0 - sqrt(f0 / g))\n",
      "        return (f0, f1, )\n",
      "\n",
      "    # Add some output to __repr__\n",
      "    def human_readable_extra(self):\n",
      "        return \"\\n\\tMulti-Objective problem\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We instantiate our problem as before, but this time we use one of the multi-objective algorithms available in PaGMO:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from PyGMO import population\n",
      "\n",
      "prob = my_mo_problem()\n",
      "algo = algorithm.sms_emoa(gen=2000)  # 2000 generations of SMS-EMOA should solve it\n",
      "pop = population(prob, 30)\n",
      "pop = algo.evolve(pop)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since in the Multi-Objective world the idea of a single 'champion' solution is not very well defined, we plot the Pareto front of the whole population, i.e., the two objectives $f_i^{(1)}$ and $f_i^{(2)}$ of each individual $i \\in 1,\\ldots,30$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "\n",
      "# Fold each objectives into vectors and print the Pareto front\n",
      "F = np.array([ind.cur_f for ind in pop]).T\n",
      "plt.scatter(F[0], F[1])\n",
      "plt.xlabel(\"$f^{(1)}$\")\n",
      "plt.ylabel(\"$f^{(2)}$\")\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "NOTE1: This problems of tutorial are implemented in PyGMO under the name **PyGMO.problem.py_example** and **PyGMO.problem.py_example_max**\n",
      "\n",
      "NOTE2: When evolve is called from an island, the process is forked and transferred to another python or ipython instance. As a consequence, when writing your **_obj_fun_impl** you cannot use stuff like matplotlib to make interactive plots and alike. If you need, during development, to have this kind of support, use the algorithm evolve method (see the optimization of the Multi-Objective problemabove)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "NOTE3: If performance is your goal, you should implement your problem in C++, and then expose it into Python."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Stochastic optimization problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let us see how to code a stochastic optimization problem, that is a problem where the objective function is stochastic.\n",
      "\n",
      "Just as before we will write a class deriving from **PyGMO.problem.base_stochastic** and reimplement some of its \u2018virtual\u2019 methods. With respect to **PyGMO.problem.base**, this base class has a field called **seed** which can be used to control the pseudo random numbers and that is changed by the algorithms that are compatible with stochastic optimization problems."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from PyGMO.problem import base_stochastic\n",
      "\n",
      "class my_problem_stochastic(base_stochastic):\n",
      "    \"\"\"\n",
      "    Noisy De Jong (sphere) function implemented purely in Python.\n",
      "\n",
      "    USAGE: my_problem_stochastic(dim = 10, seed=0)\n",
      "    \n",
      "    * dim problem dimension\n",
      "    * seed initial random seed\n",
      "    \"\"\"\n",
      "    \n",
      "    def __init__(self, dim=10, seed=0):\n",
      "        #As before, we call the base constructor, but this time with random number generator seed\n",
      "        super(my_problem_stochastic, self).__init__(dim, seed)\n",
      "        self.set_bounds(-5.12, 5.12)\n",
      "\n",
      "    def _objfun_impl(self, x):\n",
      "        from random import random as drng\n",
      "        from random import seed\n",
      "\n",
      "        #We initialize the random number generator using the\n",
      "        #data member seed (in base_stochastic). This will be changed by suitable\n",
      "        #algorithms when a stochastic problem is used. The mod operation avoids overflows\n",
      "        seed(self.seed)\n",
      "\n",
      "        #We write the objfun using the same pseudorandonm sequence\n",
      "        #as long as self.seed is unchanged.\n",
      "        f = 0.0\n",
      "        for i in range(self.dimension):\n",
      "            noise = (2 * drng() - 1) / 10\n",
      "            f = f + (x[i] + noise) * (x[i] + noise)\n",
      "        return (f, )\n",
      "    \n",
      "    def human_readable_extra(self):\n",
      "             return \"Seed: \" + str(self.seed)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us now try to solve our stochastic sphere problem using Particle Swarm Optimization algorithm"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from PyGMO import algorithm, island\n",
      "\n",
      "prob = my_problem_stochastic(dim=10, seed=123456)\n",
      "algo = algorithm.pso_gen(gen=1)\n",
      "print(\"Seed before island initialization: {}\".format(prob.seed))\n",
      "isl = island(algo, prob, 20)\n",
      "print(\"Seed after island initialization: {}\".format(prob.seed))\n",
      "for i in range(30):\n",
      "    isl.evolve(1)\n",
      "    print(\"Best F:{}, seed:{}\".format(isl.population.champion.f, isl.problem.seed))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As you can see, the current seed automatically changes during the algorithm run."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Adding a new optimization algorithm"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this section we will learn how to code a simple optimization algorithm. We will write the algorithm so that it manage multi-objective, mixed_int, constrained optimization as this will allow us to explain all the basic PyGMO functioning. Clearly our algorithm will not be very good, but a random search is always useful to benchmark against :)\n",
      "\n",
      "In a nutshell ... we will write a class deriving from PyGMO.algorithm.base and reimplement some of its \u2018virtual\u2019 methods, the main one being evolve!!!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from PyGMO.problem import base\n",
      "\n",
      "\n",
      "class my_algorithm(base):\n",
      "    \"\"\"\n",
      "    Monte-Carlo (random sampling) algorithm implemented purely in Python.\n",
      "    \"\"\"\n",
      "\n",
      "    def __init__(self, iter=10):\n",
      "        \"\"\"\n",
      "        Constructs a Monte-Carlo (random sampling) algorithm\n",
      "\n",
      "        USAGE: algorithm.my_algorithm(iter=10)\n",
      "\n",
      "        NOTE: At the end of each iteration, the randomly generated\n",
      "                point substitutes the worst individual in the population if better\n",
      "\n",
      "        * iter: number of random samples\n",
      "        \"\"\"\n",
      "\n",
      "        #We start calling the base constructor\n",
      "        super(my_algorithm, self).__init__()\n",
      "        #We then define the algorithm 'private' data members\n",
      "        self.__iter = iter\n",
      "\n",
      "    #This is the 'juice' of the algorithm, the method where the actual optimzation is coded.\n",
      "    def evolve(self, pop):\n",
      "        #If the population is empty (i.e. no individuals) nothing happens\n",
      "        if len(pop) == 0:\n",
      "            return pop\n",
      "\n",
      "        #Here we rename some variables, in particular the problem\n",
      "        prob = pop.problem\n",
      "\n",
      "        #Its dimensions (total and continuous)\n",
      "        dim, cont_dim = prob.dimension, prob.dimension - prob.i_dimension\n",
      "\n",
      "        #And the lower/upper bounds for the chromosome\n",
      "        lb, ub = prob.lb, prob.ub\n",
      "        import random\n",
      "\n",
      "        #The algorithm now starts manipulating the population\n",
      "        for _ in range(self.__iter):\n",
      "            # We create a random vector within the bounds.\n",
      "            # First the continuous part:\n",
      "            tmp_cont = [random.uniform(lb[i], ub[i]) for i in range(cont_dim)]\n",
      "\n",
      "            # Then the integer part:\n",
      "            tmp_int = [float(random.randint(lb[i], ub[i])) for i in range(cont_dim, dim)]\n",
      "\n",
      "            # Assemble them into one decision vector:\n",
      "            tmp_x = tmp_cont + tmp_int\n",
      "\n",
      "            # Push back in the population:\n",
      "            pop.push_back(tmp_x)\n",
      "\n",
      "            # Remove the worst individual:\n",
      "            pop.erase(pop.get_worst_idx())\n",
      "\n",
      "        # return the 'evolved' population\n",
      "        return pop\n",
      "\n",
      "    def get_name(self):\n",
      "        return \"Monte Carlo Algorithm (Python)\"\n",
      "\n",
      "    def human_readable_extra(self):\n",
      "        return \"n_iter=\" + str(self.__n_iter)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The above code contains a lot of interesting points worth to be discussed. So, we start\n",
      "\n",
      "* In PyGMO the decision vector (chromosome) is represented as an n-tuple. Its dimension and structure depends on the problem. Its dimension will be problem.dimension, the first prob.dimension - prob.i_dimension components will be continuous, the remaining problem.i_dimension will instead be integers.\n",
      "* When a chromosome is pushed back in a population, the domination count and the domination list (data members) are automatically updated\n",
      "* The get_worst_idx method sort the population with respect to the domination count, then domination list size (in inverse order)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}